#include "mine.h"

void fft_forward (Array<double,1> f, Array<double,1>& Fr, Array<double,1>& Fi)
{
    fftw_complex ff[J], FF[J];
    
    //Load data
    for (int j = 0; j < J; j++) {
	 ff[j][0] = f(j);
	 ff[j][1] = 0.0;
    }
    
    // Call fftw routine
    fftw_plan p = fftw_plan_dft_1d (J, ff, FF, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute (p);
    fftw_destroy_plan (p);

    //Unload data
    for (int j = 0; j < J; j++) {
	Fr(j) = FF[j][0];
	Fi(j) = FF[j][1];
    }

    // Normalize data
    Fr /= double (J);
    Fi /= double (J);
}

void fft_backward (Array<double,1> Fr, Array<double,1> Fi, Array<double,1>& f )
{
    fftw_complex ff[J], FF[J];
    
    // Load data
    for (int j = 0; j < J; j++) {
	FF[j][0] = Fr(j);
	FF[j][1] = Fi(j);
    }

    // Call fftw routine
    fftw_plan p = fftw_plan_dft_1d (J, FF, ff, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute (p);
    fftw_destroy_plan(p);

    // Unload data
    for (int j = 0; j < J; j++)
	f(j) = ff[j][0];
}

void Poisson1D (Array<double,1>& u, Array<double,1> v, double kappa)
{
    // Declare local array
    Array<double,1> Vr(J), Vi(J), Ur(J), Ui(J);
    
    //Fourier transform source term
    fft_forward (v, Vr, Vi);

    // Calculate Fourier transform of u
    Ur(0) = Ui(0) =0.0;
    for (int j = 1; j <= J/2; j++) {
	Ur(j) = -Vr(j) / double (j * j) / kappa / kappa;
	Ui(j) = -Vi(j) / double (j * j) / kappa / kappa;
    }
    for(int j = J/2; j < J; j++) {
	Ur(j) = Ur(J-j);
	Ui(j) = -Ui(J-j);
    }
    
    // Inverse Fourier transform to obtain u
    fft_backward(Ur, Ui, u);
}
